library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.3     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.0
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.0
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(glue)
library(SingleCellExperiment)
## Loading required package: SummarizedExperiment
## Loading required package: MatrixGenerics
## Loading required package: matrixStats
## 
## Attaching package: 'matrixStats'
## 
## The following object is masked from 'package:dplyr':
## 
##     count
## 
## 
## Attaching package: 'MatrixGenerics'
## 
## The following objects are masked from 'package:matrixStats':
## 
##     colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse,
##     colCounts, colCummaxs, colCummins, colCumprods, colCumsums,
##     colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs,
##     colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats,
##     colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds,
##     colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads,
##     colWeightedMeans, colWeightedMedians, colWeightedSds,
##     colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet,
##     rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods,
##     rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps,
##     rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,
##     rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks,
##     rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars,
##     rowWeightedMads, rowWeightedMeans, rowWeightedMedians,
##     rowWeightedSds, rowWeightedVars
## 
## Loading required package: GenomicRanges
## Loading required package: stats4
## Loading required package: BiocGenerics
## 
## Attaching package: 'BiocGenerics'
## 
## The following objects are masked from 'package:lubridate':
## 
##     intersect, setdiff, union
## 
## The following objects are masked from 'package:dplyr':
## 
##     combine, intersect, setdiff, union
## 
## The following objects are masked from 'package:stats':
## 
##     IQR, mad, sd, var, xtabs
## 
## The following objects are masked from 'package:base':
## 
##     anyDuplicated, aperm, append, as.data.frame, basename, cbind,
##     colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find,
##     get, grep, grepl, intersect, is.unsorted, lapply, Map, mapply,
##     match, mget, order, paste, pmax, pmax.int, pmin, pmin.int,
##     Position, rank, rbind, Reduce, rownames, sapply, setdiff, sort,
##     table, tapply, union, unique, unsplit, which.max, which.min
## 
## Loading required package: S4Vectors
## 
## Attaching package: 'S4Vectors'
## 
## The following objects are masked from 'package:lubridate':
## 
##     second, second<-
## 
## The following objects are masked from 'package:dplyr':
## 
##     first, rename
## 
## The following object is masked from 'package:tidyr':
## 
##     expand
## 
## The following object is masked from 'package:utils':
## 
##     findMatches
## 
## The following objects are masked from 'package:base':
## 
##     expand.grid, I, unname
## 
## Loading required package: IRanges
## 
## Attaching package: 'IRanges'
## 
## The following object is masked from 'package:glue':
## 
##     trim
## 
## The following object is masked from 'package:lubridate':
## 
##     %within%
## 
## The following objects are masked from 'package:dplyr':
## 
##     collapse, desc, slice
## 
## The following object is masked from 'package:purrr':
## 
##     reduce
## 
## Loading required package: GenomeInfoDb
## Loading required package: Biobase
## Welcome to Bioconductor
## 
##     Vignettes contain introductory material; view with
##     'browseVignettes()'. To cite Bioconductor, see
##     'citation("Biobase")', and for packages 'citation("pkgname")'.
## 
## 
## Attaching package: 'Biobase'
## 
## The following object is masked from 'package:MatrixGenerics':
## 
##     rowMedians
## 
## The following objects are masked from 'package:matrixStats':
## 
##     anyMissing, rowMedians
library(lemur)
## 
## Attaching package: 'lemur'
## 
## The following object is masked from 'package:dplyr':
## 
##     vars
## 
## The following object is masked from 'package:ggplot2':
## 
##     vars
source("util.R")
fit <- qs::qread("../benchmark/output/cable_spatial_plaque_data/fit_small.qs")
nei <- qs::qread("../benchmark/output/cable_spatial_plaque_data/nei_small.qs")
mouse_labels <- structure(paste0("Mouse ", 1:4), names = paste0("mouse_", 1:4))
lemur.utils::count_labels_per_neighborhood(nei, labels = vars(class %in% c("Glutamatergic", "GABAergic")), 
                                           fit = fit, cell_names = seq_len(ncol(fit))) %>%
  mutate(frac = counts / total_counts)
## # A tibble: 2 × 5
##   name  label counts total_counts     frac
##   <chr> <lgl>  <dbl>        <int>    <dbl>
## 1 Jun   TRUE    6183        30687 0.201   
## 2 Jun   FALSE      3        10368 0.000289
spatial_plaq_dens_pl <- as_tibble(colData(fit)) %>%
  ggplot(aes(x = x, y = y)) +
    ggrastr::rasterize(geom_point(aes(color = plaque_density), size = 0.05, stroke = 0), dpi = 600) +
    scale_color_gradient(low = "lightgrey", high = "darkorange", limits = c(0, 1), breaks = c(0, 0.5, 1)) +
    facet_wrap(vars(sample), labeller = as_labeller(mouse_labels)) +
    small_axis(label = "spatial coord.", fontsize = font_size_small) +
    labs(color = "Plaque density", subtitle = "Spatial structure") +
    guides(color = guide_colorbar(barheight = unit(1, "mm"))) +
    theme(legend.position = "top", strip.text = element_text(size = font_size_tiny, margin = margin(b = -1, unit = "mm")),
          plot.subtitle = element_text(size = 7))

spatial_plaq_dens_pl

set.seed(1)
umap <- uwot::umap(t(fit$embedding))
as_tibble(colData(fit)) %>%
  mutate(umap = umap) %>%
  ggplot(aes(x = umap[,1], y = umap[,2])) +
    geom_point(aes(color = sample), size = 0.3, stroke = 0) +
    guides(color = guide_legend(override.aes = list(size = 1)))

as_tibble(colData(fit)) %>%
  mutate(subclass = fct_lump(subclass, n = 5)) %>%
  mutate(umap = umap) %>%
  ggplot(aes(x = umap[,1], y = umap[,2])) +
    geom_point(aes(color = subclass == "DG"), size = 0.3, stroke = 0) +
    guides(color = guide_legend(override.aes = list(size = 1)))

as_tibble(colData(fit)) %>%
  mutate(umap = umap) %>%
  ggplot(aes(x = umap[,1], y = umap[,2])) +
    geom_point(aes(color = plaque_cluster), size = 0.3, stroke = 0) +
    guides(color = guide_legend(override.aes = list(size = 1)))

umap_plt <- as_tibble(colData(fit)) %>%
  mutate(umap = umap) %>%
  ggplot(aes(x = umap[,2], y = umap[,1])) +
    ggrastr::rasterize(geom_point(aes(color = plaque_density), size = 0.1, stroke = 0, show.legend = FALSE), dpi = 600) +
    scale_color_gradient(low = "lightgrey", high = "darkorange", limits = c(0, 1), breaks = c(0, 0.5, 1))  +
    small_axis(label = "UMAP", fontsize = font_size_small) +
    labs(subtitle = "Latent structure") +
    theme(plot.subtitle = element_text(size = 7))

umap_plt

genes_of_interest <- nei$name[1]
master_data <- as_tibble(colData(fit)) %>%
  mutate(inside = row_number() %in% nei$neighborhood[[1]]) %>%
  mutate(umap = umap) %>%
  mutate(de = assay(fit, "DE")[genes_of_interest[1], ])
make_spatial_plot <- function(data, color_by = NULL){
  data %>%
    filter(sample %in% c("mouse_1", "mouse_4")) %>%
    ggplot(aes(x = x, y = y)) +
      ggrastr::rasterize(geom_point(aes(color = {{color_by}}), size = 0.08, stroke = 0), dpi = 600) +
      facet_wrap(vars(sample), labeller = as_labeller(mouse_labels), nrow = 2) +
      small_axis(label = "spatial coord.", fontsize = font_size_small) +
      theme(strip.text = element_text(size = font_size_tiny, margin = margin(b = -1, unit = "mm")),
            legend.position = "top",
            legend.title = element_text(size = font_size),
            legend.text = element_text(size = font_size_small),
            plot.subtitle = element_text(size = 7)) 
}

make_spatial_plot(master_data, color_by = plaque_density) +
  scale_color_gradient(low = "lightgrey", high = "darkorange", limits = c(0, 1), breaks = c(0, 0.5, 1)) +
  guides(color = guide_colorbar(barheight = unit(1, "mm"), label.vjust = 3)) +
  labs(color = "")

make_umap_plot <- function(data, color_by = NULL){
  ggplot(data, aes(x = umap[,2], y = umap[,1])) +
    ggrastr::rasterize(geom_point(aes(color = {{color_by}}), size = 0.08, stroke = 0), dpi = 600) +
    small_axis(label = "UMAP", fontsize = font_size_small) +
    theme(strip.text = element_text(size = font_size_tiny, margin = margin(b = -1, unit = "mm")),
          plot.subtitle = element_text(size = 7))
}
table(fit$colData$class, fit$colData$subclass)
##                
##                 Astro CA1-ProS CA2-IG-FC   CA3    CR CT SUB    DG  Endo
##   GABAergic         0        0         0     0     0      0     0     0
##   Glutamatergic     0     9618       169  1290    53      0 13852     0
##   Non-Neuronal   3160        0         0     0     0      0     0   180
##                
##                 L2 IT ENTl L2/3 IT CTX L2/3 IT ENTl L2/3 IT PPP L2/3 IT RHP
##   GABAergic              0           0            0           0           0
##   Glutamatergic          0           0            0           1           0
##   Non-Neuronal           0           0            0           0           0
##                
##                 L4/5 IT CTX L5 PT CTX L6 CT CTX L6 IT CTX L6b CTX Lamp5
##   GABAergic               0         0         0         0       0  3082
##   Glutamatergic           0         0         0         0       0     0
##   Non-Neuronal            0         0         0         0       0     0
##                
##                 Micro-PVM NP SUB Oligo Pvalb SMC-Peri SUB-ProS  Sncg   Sst
##   GABAergic             0      0     0   208        0        0  1088   817
##   Glutamatergic         0    110     0     0        0       66     0     0
##   Non-Neuronal        640      0  6278     0       88        0     0     0
##                
##                 Sst Chodl  VLMC   Vip
##   GABAergic             0     0   333
##   Glutamatergic         0     0     0
##   Non-Neuronal          0    22     0
subclass_name_table <- c("DG" = "Dentate gyrus", "CA1-ProS" = "CA1\n(prosubiculum)", "Oligo" = "Oligodendrocytes",
                         "Astro" = "Astrocytes", "CA3" = "CA3", "CA2-IG-FC" = "CA2", "CR" = "Cajal-Retzius",
                         "Endo" = "Endothelium", "NP SUB" = "Cortex L5/6\nnear-projecting", "SUB-ProS" = "Subiculum/prosubiculum",
                          "L2/3 IT PPP" = "L2/3 IT PPP"
                         )


full_cell_types_pl <- master_data %>%
  mutate(subclass_label = case_when(
    class == "GABAergic" ~ "GABAergic Neuron",
    subclass %in% c("Micro-PVM", "VLMC", "SMC-Peri") ~ "Immune/vascular cells",
    subclass %in% names(subclass_name_table) ~ subclass_name_table[as.character(subclass)],
    TRUE ~ paste0("MISSING: ", subclass)
  )) %>%
  filter(subclass_label != "L2/3 IT PPP") %>% # It's only one cell
  mutate(subclass_label = fct_infreq(subclass_label)) %>%
  ggplot(aes(x = umap[,2], y = umap[,1])) +
    ggrastr::rasterize(geom_point(aes(color = subclass_label), size = 0.08, stroke = 0, show.legend = FALSE), dpi = 300) +
    geom_density2d(data = master_data %>% filter(inside), breaks = 0.03, contour_var = "ndensity",
                   color = "black", linetype = "dashed", linewidth = 0.6) +
    small_axis(label = "UMAP", fontsize = font_size_small) +
    guides(color = guide_legend(override.aes = list(size = 3))) +
    facet_wrap(vars(subclass_label), ncol = 6)

full_cell_types_pl

convert_interval_to_number <- function(interval){
  str_mat <- str_match(interval, ".(\\d+\\.?\\d*),(\\d+\\.?\\d*)")[,2:3]
  num_mat <- array(as.numeric(str_mat), dim = dim(str_mat))
  rowMeans(num_mat)
}

plaque_expr_relation_plt <- as_tibble(colData(fit)) %>%
  mutate(expr = logcounts(fit)[1,]) %>%
  mutate(inside = row_number() %in% nei$neighborhood[[1]]) %>%
  mutate(subclass_label = case_when(
    class == "GABAergic" ~ "GABAergic Neuron",
    subclass %in% c("Micro-PVM", "VLMC", "SMC-Peri") ~ "Immune/vascular cells",
    subclass %in% names(subclass_name_table) ~ subclass_name_table[as.character(subclass)],
    TRUE ~ paste0("MISSING: ", subclass)
  )) %>%
  filter(subclass_label != "L2/3 IT PPP") %>% # It's only one cell
  add_count(subclass_label, inside) %>%
  mutate(subclass_label = fct_infreq(subclass_label)) %>%
  summarize(expr = mean(expr), .by = c(plaque_cluster, subclass_label, inside, sample, n)) %>%
  mutate(plaque_cluster_cont = convert_interval_to_number(plaque_cluster)) %>%
  mutate(n_label_insert = ifelse(inside, "Cells inside:\\;\\;\\;", "Cells outside: ")) %>%
  mutate(n_label = paste0(n_label_insert, prettyNum(n, big.mark = "\\\\,"))) %>%
  ggplot(aes(x = plaque_cluster_cont, y = expr)) +
    geom_text(data = . %>% filter(!inside) %>% distinct(subclass_label, n_label), aes(x = 0.02, y = Inf, label = n_label), 
              vjust = 1.1, hjust = 0, size = font_size_tiny / .pt) +
    geom_text(data = . %>% filter(inside) %>% distinct(subclass_label, n_label), aes(x = 0.02, y = Inf, label = n_label), 
              vjust = 2.8, hjust = 0, size = font_size_tiny / .pt) +
    geom_smooth(aes(color = inside), span = 1.5, se = FALSE, linewidth = 2, show.legend = FALSE) +
    geom_point(aes(color = stage(inside, after_scale = colorspace::darken(color, 0.2))), size = 0.3) +
    scale_color_manual(values = c("TRUE" = "darkblue", "FALSE" = "darkgrey")) +
    scale_x_continuous(limits = c(0, 1), breaks = c(0, 0.5, 1), labels = c("0", "0.5", "1"),
                       expand = expansion(add = 0)) +
    coord_cartesian(ylim = c(0, 0.5)) +
    facet_wrap(vars(subclass_label), scales = "fixed", ncol = 6) +
    guides(color = guide_legend(override.aes = list(size = 2))) +
    theme(panel.spacing.x = unit(4, "mm"), legend.position = "bottom") +
    labs(color = "Inside \\emph{Jun}'s\nneighborhood",
         x = "Amyloid-$\\beta$ density",
         y = "\\emph{Jun} expression (pseudobulked)")
plot_assemble(
  add_text("(A) Alzheimer's disease mouse model cell type details", x = 0.5, y = 1, fontsize = font_size, vjust = 1, fontface = "bold"),
  add_plot(full_cell_types_pl, x = 8,  y = 3, width = 162, height = 100),

  add_text("(B) Relation of Amyloid-$\\beta$ density and \\emph{Jun} expression per cell type", x = 0.5, y = 104, fontsize = font_size, vjust = 1, fontface = "bold"),
  add_plot(plaque_expr_relation_plt, x = 0,  y = 108, width = 170, height = 66),
  
  
  width = 170, height = 176, units = "mm", show_grid_lines = FALSE,
  latex_support = TRUE, filename = "../plots/suppl_alzheimer_cell_types.pdf"
)
## Using TikZ metrics dictionary at:
##  alzheimer_plaques_analysis-tikzDictionary
## gg[gg1]
## gg[gg2]
## gg[gg3]
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : span too small.  fewer data values than degrees of freedom.
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : pseudoinverse used at 0.0455
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : neighborhood radius 1.1078
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : reciprocal condition number 0
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : There are other near singularities as well. 1.2272
## gg[gg4]
##  [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
amb_sp <- make_spatial_plot(master_data, color_by = plaque_density) +
  scale_color_gradient(low = "lightgrey", high = "darkorange", limits = c(0, 1), 
                       breaks = c(0, 0.5, 1), labels = c("0", "0.5", "1")) +
  guides(color = guide_colorbar(barheight = unit(1, "mm"), label.vjust = 3)) +
  labs(color = "")
amb_umap <- make_umap_plot(master_data, color_by = plaque_density) +
  scale_color_gradient(low = "lightgrey", high = "darkorange", limits = c(0, 1), breaks = c(0, 0.5, 1)) +
  guides(color = "none")


ct_master_data <- master_data %>%
  mutate(is_neuron = ifelse(class %in% c("Glutamatergic", "GABAergic"), "neuron", "other"))
ct_sp <- make_spatial_plot(ct_master_data, color_by = is_neuron) +
  scale_color_manual(values = c("neuron" = "lightblue", "other" = "lightgrey"), labels = c("neuron"="Neuron", "other"="other cell type")) +
  guides(color = guide_legend(override.aes = list(size = 3))) +
  labs(color = "") +
  theme(legend.text = element_text(size = font_size))
ct_umap <- make_umap_plot(ct_master_data, color_by = is_neuron) +
  scale_color_manual(values = c("neuron" = "lightblue", "other" = "lightgrey")) +
  guides(color = "none")

abs_max <- max(abs(quantile(assay(fit, "DE")["Jun", ], prob = c(0.05, 0.95))))
de_sp <- make_spatial_plot(master_data, color_by = de) +
  guides(color = guide_colorbar(barheight = unit(1, "mm"), label.vjust = 3)) +
  scale_color_de_gradient(abs_max = 0.2, breaks = c(-0.2, 0, 0.2), labels = c("-0.2", "0", "0.2")) +
  labs(color = "")
de_umap <- make_umap_plot(master_data, color_by = de) +
  scale_color_de_gradient(abs_max = abs_max) +
    guides(color = "none")

nei_master_data <- master_data %>%
  mutate(inside = ifelse(inside, "in", "out"))
nei_sp <- make_spatial_plot(nei_master_data, color_by = inside) +
  scale_color_manual(values = c("in" = "darkblue", "out" = "lightgrey"), labels = c("in"="Inside \\emph{Jun} neighborhood", "out"="outside")) +
  guides(color = guide_legend(override.aes = list(size = 3))) +
  labs(color = "") +
  theme(legend.text = element_text(size = font_size))
nei_umap <- make_umap_plot(nei_master_data, color_by = inside) +
  scale_color_manual(values = c("in" = "darkblue", "out" = "lightgrey")) +
  guides(color = "none")
genes_of_interest <- nei$name[1]
mask <- matrix(NA, nrow = length(genes_of_interest), ncol = ncol(fit), 
               dimnames = list(genes_of_interest, colnames(fit)))
mask2 <- matrix(NA, nrow = length(genes_of_interest), ncol = ncol(fit), 
               dimnames = list(genes_of_interest, colnames(fit)))
mask3 <- matrix(NA, nrow = length(genes_of_interest), ncol = ncol(fit), 
               dimnames = list(genes_of_interest, colnames(fit)))


for(id in genes_of_interest){
  mask[id, filter(nei, name == id)$neighborhood[[1]]] <- 1
  
  mask2[id, ! fit$colData$class %in% c("Glutamatergic", "GABAergic")] <- 1
  mask2[id, filter(nei, name == id)$neighborhood[[1]]] <- NA
  
  mask3[id, fit$colData$class %in% c("Glutamatergic", "GABAergic")] <- 1
  mask3[id, filter(nei, name == id)$neighborhood[[1]]] <- NA
}

psce2 <- glmGamPoi::pseudobulk(SingleCellExperiment(list(inside = as.matrix(logcounts(fit)[genes_of_interest,,drop=FALSE] * mask),
                                                        outside = as.matrix(logcounts(fit)[genes_of_interest,,drop=FALSE] * mask2),
                                                        outside_neurons = as.matrix(logcounts(fit)[genes_of_interest,,drop=FALSE] * mask3))),
                      group_by = vars(sample, plaque_cluster), n = n(),
                      aggregation_functions = list(.default = \(...) matrixStats::rowMeans2(..., na.rm = TRUE)),
                      col_data = as.data.frame(colData(fit)))
## Aggregating assay 'inside' using 'custom function'.
## Aggregating assay 'outside' using 'custom function'.
## Aggregating assay 'outside_neurons' using 'custom function'.
comparison_data <- as_tibble(colData(psce2)) %>%
  mutate(expr_inside = as_tibble(t(assay(psce2, "inside"))),
         expr_outside = as_tibble(t(assay(psce2, "outside"))),
         expr_outsideNeuron = as_tibble(t(assay(psce2, "outside_neurons")))) %>%
  unpack(starts_with("expr"), names_sep = "-") %>%
  pivot_longer(starts_with("expr"), names_sep = "[-_]", names_to = c(".value", "origin", "symbol")) 



expr_comparison_pl <- comparison_data %>%
  mutate(plaque_cluster_mean = convert_interval_to_number(as.character(plaque_cluster))) %>%
  mutate(origin = factor(origin, levels = c("inside", "outsideNeuron", "outside")))  %>%
  ggplot(aes(x = plaque_cluster_mean, y = expr)) +
    geom_point(aes(color = origin), size = 0.3) +
    geom_smooth(aes(color = origin), span = 1.5, se = FALSE, linewidth = 2) +
    scale_color_manual(values = c("inside" = "darkblue", "outsideNeuron" = "lightblue", "outside" = "lightgrey"),
                       labels = c("inside" = "Cells in neighborhood", "outsideNeuron" = "Other neurons", "outside" = "All other cells")) +
    scale_y_continuous(limits = c(0, 0.3), expand = expansion(add = 0)) +
    scale_x_continuous(breaks = seq(0, 1, by = 0.1), limits = c(0, 1)) +
    theme(legend.position = "right") +
    labs(color = "", y = "Expression", x = "Amyloid-$\\beta$ plaque density (binned)",
         subtitle = "\\emph{Jun} expr. vs. Amyloid-$\\beta$ plaque density")  +
    theme(plot.subtitle = element_text(size = font_size))

expr_comparison_pl
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

plot_assemble(
  add_text("(A) ", x = 0.5, y = 1, fontsize = font_size, vjust = 1, fontface = "bold"),
  add_text("Amyloid-$\\beta$ density", x = 6, y = 2, fontsize = font_size, vjust = 1, fontface = "plain"),
  add_plot(my_get_legend(amb_sp), x = 37, y = 3, height = 5, width = 1),
  add_plot(amb_sp + guides(color = "none"),   x = 2,   y = 5, width = 23, height = 45),
  add_plot(amb_umap, x = 22,  y = 5, width = 33, height = 45),

  add_text("(B) ", x = 56, y = 1, fontsize = font_size, vjust = 1, fontface = "bold"),
  add_plot(my_get_legend(ct_sp), x = 84, y = 0.5, height = 5, width = 1),
  add_plot(ct_sp + guides(color = "none"),    x = 58,  y = 5, width = 23, height = 45),
  add_plot(ct_umap,  x = 78,  y = 5, width = 33, height = 45),

  add_text("(C) ", x = 112, y = 1, fontsize = font_size, vjust = 1, fontface = "bold"),
  add_text("\\emph{Jun} differential expression", x = 118, y = 2, fontsize = font_size, vjust = 1, fontface = "plain"),
  add_plot(my_get_legend(de_sp), x = 159, y = 3, height = 5, width = 1),
  add_plot(de_sp + guides(color = "none"),    x = 114, y = 5, width = 23, height = 45),
  add_plot(de_umap,  x = 134, y = 5, width = 33, height = 45),

  add_text("(D) ", x = 0.5, y = 50, fontsize = font_size, vjust = 1, fontface = "bold"),
  add_plot(nei_sp,   x = 2,   y = 50, width = 23, height = 50),
  add_plot(nei_umap, x = 22,  y = 55, width = 33, height = 45),
  
  add_text("(E) ", x = 65, y = 50, fontsize = font_size, vjust = 1, fontface = "bold"),
  add_plot(expr_comparison_pl, x = 66, y = 50, width = 100, height = 50),
  
  width = 170, height = 100, units = "mm", show_grid_lines = FALSE,
  latex_support = TRUE, filename = "../plots/alzheimer_figure.pdf"
)
## gg[gg1]
## gg[gg2]
## gg[gg3]
## gg[gg4]
## gg[gg5]
## gg[gg6]
## gg[gg7]
## gg[gg8]
## gg[gg9]
## gg[gg10]
## gg[gg11]
## gg[gg12]
## gg[gg13]
## gg[gg14]
## gg[gg15]
## gg[gg16]
## gg[gg17]
## gg[gg18]
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
## gg[gg19]
##  [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
master_data %>%
  mutate(label = case_when(
    inside ~ "inside",
    cell_type == "NEURON" & ! inside ~ "other neurons",
    TRUE ~ "other",
  )) %>%
  sample_frac(n = 1) %>%
  ggplot(aes(x = x, y = y)) +
    geom_point(aes(color = label), size = 0.1) +
    facet_wrap(vars(sample))